national_data$yes[national_data$department_id==1468] <- 1

# bivariate
model_bivariate <- lm(yes~condition, data=national_data)
model_bivariate_robust <- coeftest(model_bivariate, vcov.=vcovHC(model_bivariate, type="HC1")) #[1:2,]

#cluster
model_cluster <- lm(yes~condition+factor(match_id), data=national_data)

model_cluster_robust <- coeftest(model_cluster, vcov.=rddtools::vcovCluster(model_cluster, clusterVar=national_data$match_id), cluster = national_data$match_id) #[1:2,]


#cluster and rank
model_cluster_rank <- lm(yes~condition*rank + factor(match_id), data=national_data)

model_cluster_rank_robust <-coeftest(model_cluster_rank, vcov.=rddtools::vcovCluster(model_cluster_rank, clusterVar=national_data$match_id), cluster = national_data$match_id)

#cluster and rank tercile
model_cluster_rank_tercile <- lm(yes~condition*ranktercile + factor(match_id), data=national_data)

model_cluster_rank_tercile_robust <-coeftest(model_cluster_rank_tercile, vcov.=rddtools::vcovCluster(model_cluster_rank_tercile, clusterVar=national_data$match_id), cluster = national_data$match_id) #[c(1:4, 1251:1252),]

gaze.coeft <- function(x, col="Std. Error"){
	stopifnot(is.list(x))
	out <- lapply(x, function(y){
		y[ , col]
	})
	return(out)
}

stargazer(model_bivariate,model_cluster, model_cluster_rank, model_cluster_rank_tercile,
		  type = "latex", 
		  se = gaze.coeft(list(model_bivariate_robust, model_cluster_robust, model_cluster_rank_robust,model_cluster_rank_tercile_robust)), 
		  omit = "match_id",
		  add.lines = list(c("Matched Pair FE", "No", "Yes", "Yes", "Yes")),
		  intercept.bottom = FALSE, 
		  keep.stat = c("rsq", "n"), 
		  digits=2, 
		  star.cutoffs = c(0.05, 0.01, 0.001),
		  dep.var.labels = "Agreed to Discuss Collaboration",
		  covariate.labels = c("Intercept", "Ranking Condition", "Numeric Rank","Middle Rank Tercile", "Top Rank Tercile", "Ranking Condition X Numeric Rank","Ranking Condition X Middle Rank Tercile","Ranking Condition X Top Rank Tercile"), out = "results/tablej1.tex")